Nature Genetics|从GWAS结果里挖掘候选基因方法|群体遗传专题

进行动植物性状的基因定位是功能基因组时代的不变热点。GWAS是进行性状定位的有力手段之一，也被广泛的应用于各种物种中。确定GWAS关联区间一般有两种方法，一种是peak SNP上下游200kb，这个方法简单粗暴，每个位点的区间也是固定不变。一种是依赖于LD block，由于染色体不同位置的LD不同，所以这种方法区间大小都不一致，得到的区间大小也更准确。当然，GWAS也可以和遗传图谱QTL联合分析，具体内容详见上上期的2020 PBJ|关联分析和连锁分析解析玉米籽粒大小的遗传结构|群体遗传专题。本期我们聚焦于基于GWAS研究结果筛选可靠的候选基因。

在筛选候选基因之前，我们先聊下候选基因的常规筛选方法。在DNA层面，多数情况下，功能基因非同义突变会引起表型的差异，所以通过筛选定位区间的非同义突变是一种鉴定候选基因的常见策略。但是也有一些例外的情况，比如promoter区的突变，或者内含子区的突变也有可能是潜在的causal基因。RNA层面，差异表达基因的筛选也是一种常见的策略，而且可以和非同义突变的候选基因筛选联合起来分析。由于基因变异的复杂性，causal基因的变异也有其他的类型，比如大片段的插入缺失等等。本次我们通过一篇2016年online的Nature Genetics经典文章解析一种比较新的候选基因筛选策略。同样以材料和方法、结果和靠谱er评价结构展开。

材料和方法

作者选择了176个来自日本粳稻品种。现在看来材料数量偏少。然后进行重测序，平均测序深度为5.8 X，现在看来测序深度也偏低。基于BWA比对和GATK SNP calling。SNP完整度和MAF过滤标准分别为0.75和0.05。转录起始位点的上游2Kb定义为promoter区。

通过EIGENSTRAT进行PCA分析，PHYLIP邻接法构建进化树。PLINK进行LD decay分析，R包可视化。LD的r²＞0.6作为候选区间筛选标准。LMM模型进行GWAS关联分析。关联阈值采用false-discovery-rate-adjusted P values。模拟1000次。候选区间的标记用1和-1进行编码，1和-1分别代表主要多态和稀有多态。单倍型可视化采用bin的方式进行，bin大小为10Kb。每个bin通过K-means聚类，并强行聚成2个group。其他方法详见原文。

结果

由于本文的材料方法部分并未描述具体的候选基因鉴定的方法，我们尊重原文的安排。但是为了让大家在看结果之前对本文的方法有个清晰的认识，我们将本文的候选基因筛选的核心方法放在结果描述部分之前。如下所示：

首先，将GWAS候选区间的所有SNP按照以下的规则进行分类（共五类）。

I：引起氨基酸编码变化和可变剪切位点（位于内含子start或者end处的GT或者AG）的显著关联SNP；

II：位于基因5’端的显著关联SNP（比如起始密码子ATG上游2Kb内的SNP）；

III：不满足以上I和II标准的显著关联SNP（比如位于编码区，但是并未引起氨基酸变化的内含子区或者3’端非编码区）；

IV：编码区（coding regions）外的显著关联SNP；

V：不显著关联的SNP。

下一步就基于分类的结果，重点提取group I的SNP进行单倍型分析，结合不同单倍型的表型判断潜在的causal单倍型，然后进行基因组注释、GO和KEGG注释和/或同源基因比对筛选候选基因。如果group I中筛选不到潜在的候选基因则考虑其他的group内的SNP。比如group IV，有些基因间区的lncRNA也会对表型有调控作用，依此类推。

以上的方法也可以反着进行，比如先进行group I的SNP筛选，并根据各种注释筛选到潜在的候选基因，然后进行单倍型分析。 

首先本文的材料，作者花了一段的篇幅说明了这176个品种的多样性较高，并说了选择这176个品种的理由，这些内容的潜台词似乎就是要告诉读者：我这是指哪儿打哪儿。

本文有偏方法论的内容，所以为了使方法更robust，作者的表型重点选择了抽穗期，这个表型被克隆的基因有很多，方便作者作为内参/阳性对照使用。

针对抽穗期的定位，在不同染色体上获得了多个显著关联位点，其中第3号，6号和7号染色体上的loci分别有已经被克隆的基因。针对7号染色体上的loci，该位点附近有已经克隆的基因Hd2。LD block分析将该位点锁定在346Kb区间，该区间有124个多态性位点，其中只有1个位点位于group I，该位点的基因就是已经克隆的基因Hd2。该变异包含了3种单倍型，其中C型单倍型和Hd2一致。通过Hd2基因的定位说了该方法的准确性。

对新基因的挖掘，以1号染色体上的位点为例，通过LD block将该位点锁定在346Kb的区间，其中位于group I的变异中有7个基因，其中一个基因LOC_Os01g62780 是拟南芥中的HESO1的同源基因，该基因突变在拟南芥中展示出延迟开花的特征。单倍型分析发现该基因有两种单倍型（A和B），其中B单倍型的SNP位于I组，B单倍型的抽穗期比A型抽穗期延迟。克隆A和B型基因并转入A型材料日本晴中，发现转A型基因材料中抽穗期没有明显变化，而转B型基因材料中表现出了明显的延迟抽穗。因此证明了LOC_Os01g62780就是抽穗期的候选基因。

用类似的套路对11号染色体上的新基因进行了鉴定。

总结

这种方法看上去对鉴定候选非常有效，但是该方法可行的遗传学内核仍然在于连锁不平衡。本文判断连锁不平衡的依据是LD 的r2＞0.6。如果给一个更严格的r2阈值，我们的候选区间会更小。这里靠谱er给个不靠谱的私货：理论上单倍型的大小就是LD的极限。这也就是为什么进行LD分析又能继续进行单倍型分析筛候选基因的原因所在（如果有误，大神请轻点拍砖），但是在确定候选区间时，需要在定位的精度和准度之间寻找一个平衡，单纯追求精度会导致causal基因被漏掉的概率大大增加。因此，实际应用中，会对LD block阈值大小做个权衡（trade off）。第二点，如果针对一个参考基因组不那么完善的物种，要获得单倍型信息就需要进行PCR一代测序了（或者基于三代测序平台的重测序）。第三点，等位基因的异质性可能会导致定位的结果产生偏离。这一点作者也在原文中进行了讨论，因为等位基因过多会导致统计学上显著性丧失。这一点可以通过增加样本量或者进行全基因组gene-based 关联分析解决。

参考文献

Yano K, Yamamoto E, Aya K, et al. Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice[J]. Nature genetics, 2016, 48(8): 927.

本文作者-冯胜

冯胜，负责海南&湛江地区业务，硕士毕业后入职联川，已服务近200客户，顺利完成200+项目。

已经合作客户单位中国热带农业科学院、海南大学、海南师范大学、海南医学院、海南省农业科学院、广东海洋大学、广东医科大学（湛江校区）、海南省人民医院、海南省中医院、海口市人民医院、广东医科大学附属医院等。

入职联川的2年多时间里，在联川浓郁的学习氛围中，从一个科研小白迅速成长为技术型销售。

学习

基因组测序作为群体分子育种的常规手段之一，入职联川之前，作为一个连基因组是什么都不知道的小白，现对群体遗传进化、遗传图谱、GWAS、BSA和BSR项目从售前沟通、群体选择、基本项目方案，在联川群体大佬许语辉博士的小灶教学下，得到了全面的提升，可以做到无障碍沟通。

专业

参编行业首本群体遗传书籍-《NGS时代的BSA百科全书》，从实验设计、样本选取、测序技术、分析流程到常见问题解决方案，均有专业详细解释，研一小白看透了都可以轻松地设计出完整项目方案。

责任

责任这一个词在联川体现的淋淋尽致，20年4月底客户省基金申报，凌晨12点多电话联系沟通项目方案，微信留言技术支持，二话不说打开电脑，凌晨一点准时将项目方案交付客户（这种做法不建议哦）。

作为一名专业型销售，从售前的实验方案设计到售后的数据整理、指导云平台个性化在线作图，我会全程跟进服务。项目交给我请您放心！

相关阅读